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Abstract. The radiation (reaction, Robin) boundary condition for the diffusion equation is widely used 
in chemical and biological applications to express reactive boundaries. The underlying trajectories of the 
diffusing particles are believed to be partially absorbed and partially reflected at the reactive boundary, 
however, the relation between the reaction constant in the Robin boundary condition and the reflection 
probability is not well defined. In this paper we define the partially reflected process as a limit of the Marko- 
vian jump process generated by the Euler scheme for the underlying Ito dynamics with partial boundary 
reflection. Trajectories that cross the boundary are terminated with probability P\/ At and otherwise are 
reflected in a normal or oblique direction. We use boundary layer analysis of the corresponding master 
equation to resolve the non-uniform convergence of the probability density function of the numerical scheme 
to the solution of the Fokker-Planck equation in a half space, with the Robin constant k. The boundary 
layer equation is of the Wiener-Hopf type. We show that the Robin boundary condition is recovered if 
and only if trajectories are reflected in the co-normal direction crn, where <r is the (possibly anisotropic) 
constant diffusion matrix and n is the unit normal to the boundary. Otherwise, the density satisfies an 
oblique derivative boundary condition. The constant k is related to P by n = rP^/un, where r = l/V^ an d 
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1. Introduction. The Fokker-Planck equation (FPE) with radiation (also called reac- 
tive or Robin) boundary conditions is widely used to describe diffusion in a biological cell 
with chemical reactions on its surface [T], [2], [3], [J, [5], [6], [7], [8], [9]. The Robin bound- 
ary conditions are used in [2] , [4] , [5] , [6] as a homogenization of mixed Dirichlct- Neumann 
boundary conditions given on scattered small absorbing windows in an otherwise reflecting 
boundary. The latter may represent, e.g., ligand binding or pumping out ions at sites on the 
boundary of a biological cell and no flux through the remaining boundary. The reactive rate 
constant in the Robin boundary conditions is chosen in the homogenization process so that 
the decay rate of the survival probability is the same as that in the mixed Dirichlet- Neumann 
boundary value problem. 

The definition of the Ito stochastic dynamics 



x = a(x, t) + y^2a(x, t) w, (I- 1 ) 

on the positive axis with total or partial reflection at the origin was given first by Feller 
[10] for the one-dimensional case with a(x, t) and a{x,t) independent of i, as a limit of Ito 
processes, which are terminated when they reach the boundary or moved instantaneously 
to a point x = Pj > with probability p« . When pj — > 1 and pj — > with 

lim — = c, (1.2) 

j'^oo pj 

where c is a constant, the partially reflected process converges to a limit. The transition prob- 
ability density function (pdf) of the limit process, p(y, t \ x, s) dy — Pr{x(t) £ (y,y + dy) \ x(s) — x}, 
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is the solution of the FPE 

dp(y,t\x,s) _ d[a(y,t)p(y,t\x,s)} d 2 [<r(y,t)p(y,t\x,s)} 
dt dy dy 2 

or equivalently, 

dp(y,t\x,s) = _dJ(y,t\x,s) forall y ^ x> ^ 

where 

J(y, t 1 x, s) = a(y, t) P (y, t 1 x, s) - ^)p(y,t\x,s)] ^ 

is the flux. The initial condition is 

p(y,t | x, s) — * S(y — x) as t[s (1.5) 

and the radiation boundary condition is 

- J(0,t \x,s) = Kp(0,t\x,s), (1.6) 

where k is a constant related to the constant c and to the values of the coefficients at the 
boundary. The no flux and Dirichlct boundary conditions are recovered if c = or c = oo, 
respectively. Feller's method does not translate into a Brownian dynamics simulation of 
the limit process, because his approximations are continuous-time Ito processes. Skorokhod 
[TT] defines the reflection process inside the boundary. Several numerical schemes have been 
proposed for simulating this process (see, e.g., [TT], [12], [13], [H]). The main issue there is 
to approximate the local time spent on the boundary. 

The definition of a diffusion process with absorbing or reflecting boundaries as limits of 
Markovian jump processes, which is the basis for all simulations, gives in the limit diffusion 
processes with well defined boundary behavior. However, the definition of a diffusion process 
with partially reflecting boundaries as a limit of Markovian jump processes gives different 
diffusions for different jump processes. This is expressed in different relations between the 
termination probability of the jump process and the boundary conditions for the Fokker- 
Planck equations (see, e.g., [H]). The process x(t) defined by equation (jl.ip with partially 
absorbing boundaries can be defined as the limit of the solutions of the Markovian jump 
processes generated by the Euler scheme 

x At (t + At) = x At (t) +a(x At {t),i) At + \/2a(x At (t),t) Aw(t,At) for t>s (1.7) 



(1.3) 



(1.4) 



XAt(s) = x (1.8) 

in the interval x > 0, for < t ~ s < T, with At = T/N, t- s = iT/N (i = 0, 1, . . . , N), 
where for each t the random variables Aw(t, At) are normally distributed and independent 
with zero mean and variance A£. The partially absorbing boundary condition for (|1.7|) has 
to be chosen so that the pdf p A t(x, t) of x A t(t) converges to the solution of (|1.3[) - (|1.6p . At a 
partially reflecting boundary for (|1.7jl . the trajectories are reflected with probability (w.p.) 
R and otherwise terminated (absorbed), once they cross the origin. We show below that 
keeping R constant (e.g., R = 1/2) as At — > leads to the convergence of the pdf p A t(x,t) 
to the solution of the FPE with an absorbing rather than the Robin boundary condition. 
Thus the reflection probability R must increase to 1 as At — > in order to yield the Robin 
condition (|1.6[) . Moreover, the reactive constant k is related to the limit 

1-R , . 

lim —= = P. 1.9) 
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The reflecting boundary condition is recovered for P = 0, while the absorbing boundary 
condition is obtained for P = oo. Motivated by these considerations, we design the following 
simple boundary behavior for the simulated trajectories that cross the boundary, identified 
by XAt(t) + a(x At (t),t)At + y/2a(x A t(t),t) Aw < 0, 

{ -(x A t(t) + a(x At (t),t)At + y/2a(x At (t),t) Aw) w.p. 1 - P^/At 
x At {t + At) = ! (1.10) 
[ terminate trajectory otherwise. 

The exiting trajectory is normally reflected w.p. 

R=l-PVAi (1.11) 

and is otherwise terminated (absorbed). The scaling of the termination probability with y/ At 
reflects the fact that the discrete unidirectional diffusion current at any point, including the 

boundary, is O (l/y/Atj (see [15], QS]). This means that the number of discrete trajectories 

hitting or crossing the boundary in any finite time interval increases as 1 /V At. Therefore, 
to keep the efflux of trajectories finite as At — > 0, the termination probability of a crossing 
trajectory, 1 — R, has to be 0(y/~A~i). The pdf p A t(x, t), however, does not converge to the 
solution p(x,t) of (|1.3[) - (|1.6[) on the boundary, as discussed in Section [21 This is due to the 
formation of a boundary layer, as is typical for diffusion approximations of Markovian jump 
processes that jump over the boundary [17] . [18] , [19j . The boundary layer equations are 
typically Wiener- Hopf integral equations. The Wiener-Hopf boundary layer equation for the 
particular case of a partially reflected Brownian motion on the positive axis (i.e., a(x, t) = 
and <r(x,t) = a in (|1.7p ) was recently solved in [5j and the relationship k = P\fa / ' ypK was 
found. 

The convergence of the pdf of an Euler scheme has been studied in [5D], [5T] for the 
higher-dimensional problem with oblique reflection. Bounds on the integral norm of the 
approximation error are given for the solution of the backward Kolmogorov equation. These, 
however, do not resolve the boundary layer of the pdf of the numerical solution. The solution 
of the forward equation for the Euler scheme converges non-uniformly to the solution of the 
Fokker-Planck equation due to the appearance of a boundary layer in the first order spatial 
derivative. This distorts the boundary flux and gives incorrect boundary conditions. A 
boundary layer expansion is needed to capture the boundary phenomena. 

The derivation of the radiation condition has a long history. Collins and Kimball 
[22] (see also [23]) derived the radiation boundary condition (|1.6p for the limit p(x,t) = 
\im. At —>oPAt(x,t) from an underlying discrete random walk model on a semi-infinite one- 
dimensional lattice with partial absorbtion at the endpoint. Their model assumes constant 
diffusion coefficient and vanishing drift, for which they find the reactive constant in terms 
of the absorbtion probability and the diffusion coefficient. Previous simulation schemes that 
recover the Robin boundary condition [T], [53], [25], [2B], [57] make use of the explicit so- 
lution to the half space FPE with linear drift term and constant diffusion coefficient with 
a Robin condition. In [28, and references therein] the specular reflection method near a 
reflecting boundary has been shown to be superior to other methods, such as rejection, 
multiple rejection and interruption. 

An apparent paradox arises when using (jl.7p and other schemes: while the pdf p A t{y, t\x,s) 
of the solution of (fT77|) . fTS"]) . (fTTTDj) . (fTTTT]) converges to the solution of the FPE (fOfl and 
the initial condition (|1.5[) . each approximant p At (y,t\x, s) does not satisfy the boundary 
condition (|1.6p , not even approximately, that is, the error does not decay as At — > 0. For a 
general diffusion coefficient and drift term, the boundary condition is not satisfied even for 
the case of a reflecting boundary condition. This problem plagues other schemes as well. 
The apparent paradox is due to the non- uniform convergence of p A t(y, t\x,s) to the solution 
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p(y, t\x,s) of the Fokker-Planck equation, caused by a boundary layer in PAt(y, t 1 x , s), as is 
typical of boundary behavior of diffusion approximations to Mar kovian jump processes. The 
limit p(y,t \x,s), however, satisfies the boundary condition (|f .6[) for some k. Our analysis 
can be extended to other schemes in a straightforward way. ft is well known that the Euler 
scheme produces an 0(a/A?) error in estimating the mean first passage time to reach an ab- 
sorbing boundary. There are several recipes to reduce the discretization error to O(At) [29] . 
|30j , |31j , [32j , [33] . Another manifestation of the boundary layer is that the approximation 
error of the pdf near absorbing or reflecting boundaries is 0{\fKt), and methods, including 
PQ, [33] reduce this error to O(At). Thus, we expect the formation of a boundary layer of 
size 0(%/A£) for the Euler scheme (|1.7p with the boundary behavior (|1.10|) . 

This paper is concerned with the convergence of the partially reflecting Markovian jump 
process generated by (|1.7p . (| 1 . lOf) in one and higher dimensions. We show that this scheme, 
with the additional requirement that the pdf converges to the solution of the FPE with a 
given Robin boundary condition, defines a unique diffusion process with partial reflection 
at the boundary. This definition is then generalized to higher dimensions. In contrast 
to the Collins and Kimball [25] discrete scheme, this definition is not restricted to lattice 
points and the drift and diffusion coefficients may vary. The advantage of the current 
suggested design (| 1 . 10[) is its simplicity, which is both easily and efficiently implemented 
and amenable to analysis. There is no need to make any assumptions on the structure of 
the diffusion coefficient or the drift. From the theoretical point of view, it serves as a physical 
interpretation for the behavior of diffusive trajectories near a reactive boundary. 

Our main result in the one-dimensional case is the relation between the reactive "con- 
stant" K,(t) and the absorbtion parameter P for the dynamics (|l.f I) on the positive axis with 
drift and with a variable diffusion coefficient, 

K(t)=rPy/a(0,t), r=-=, (1.12) 

V 7r 

The relation (|1.12[) is new for diffusion with variable coefficients. The value r = 1/y/ir is 
different than values obtained for other schemes, e.g., than the value r = 1/V2, predicted by 
the discrete random walk theory of radiation boundaries [22] . Values of r for other schemes 
are given in [8]. We show the effect of using (j!.12|) in numerical simulations. 

The scheme (11.101) is generalized to diffusion with drift and anisotropic constant diffusion 
matrix <r(t) in the half space, x\ > 0, with partial oblique reflection. We show that the 
Robin boundary condition is recovered if and only if trajectories are reflected in the direction 
of the unit vector 

an 

v = - -, 1.13 

\\trn\\ 

where n is the unit normal to the boundary. The radiation parameter k{x, t) in the d- 
dimcnsional Robin boundary condition and the absorbtion parameter P{x) are related by 



K (x,t) = rP(x)y/a n {t), xi =0, (1.14) 

with r given in (|1.12[) and a n (t) — n T cr(t)n. The relation (|1.14[) is new for higher- 
dimensional diffusion in a half space with drift and anisotropic diffusion matrix. 

In the most common case of constant isotropic diffusion our result extends to domains 
with curved boundaries. This is due to the fact that a smooth local mapping of the domain 
to a half space with an orthogonal system of coordinates preserves the constant isotropic 
diffusion matrix, though the drift changes according to Ito's formula. In this case the vector 
v coincides with the normal n. 
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2. Boundary layer analysis in one dimension. The aim of the boundary layer 
analysis below is to examine the convergence of the pdf PAtiVi t\x, s) of the solution XAt(t) 
of (|1.7p . (|1.8|) to the solution p(y, t\x,s) of ()1.3[) - (|1.6p . and to find the relation between the 
parameter P of (jl.lOp and the reactive constant k in (II. 6|) . Using abbreviated notation, 
the pdf PAt(y, 1 1 x, s) — pAt(y, t) satisfies the forward Kolmogorov equation [15] , [16] , [17], 



P At(y,t + At) 



PAt{x,t) 

o ^/Ana{x,t)At 



exp 



(y — x — a(x, t)Aty 
Aa(x,t)At 



+ 



(1 - PVAi) exp 



(y + x + a(x,t)Aty 



4a(x,t)At 

For P = the pdf pAt{y,t) satisfies the boundary condition 

dp At (0,t) 



dx. 



dy 



0. 



(2.1) 



(2.2) 



which is obtained by differentiation of (|2.1|) with respect to y at y = 0. If P ^ 0, we obtain 

d P At(0,t + At) PAt (0,t)P 



dy 

which holds also in the limit At - 
matters, indeed, 



rO(VAt), 

0. However, the order of the limits At 



lixn lim^'^lim lim *^M. 
At^O yio dy yio At^o dy 



(2.3) 
and y i 

(2.4) 



The limit of (12. 3|) is not the boundary condition that the limit function p{y 1 1) = limAt^oPAt(j/, t) 
(for y > 0) satisfies. To find the boundary condition of p(y, t), in either case P = or P ^ 0, 
we show below that p(y, t) satisfies the FPE (|1.3|) and the initial condition (|1.5j) for all y > 0. 
Since for P = the simulation preserves probability (the population of trajectories), 



d f°° 

0= — / p(x,t)dx = — 
dtJ Q Py ' 



d[<T(0,t)p(0,t)} 

dy 



■a{0,t)p(0,t) = J(0,t). 



(2.5) 



Equation (|2.5p is the no-flux boundary condition. The discrepancy between (|2.5[) and 
is due to the nonuniform convergence of PAt(y,t) to its limit p(y,t) in the interval. There 
is a boundary layer of width 0(\/Ai), in which the boundary condition (|2.2[) for pAt(y,t) 
changes into the boundary condition (|2.5[) that p(y, t) satisfies. To analyze the discrepancy 
between (12.21) and (|2.5p . we introduce the local variable y = -q^/At and the boundary layer 
solution 

p BL (v,t) =PAt(v^Al,t)- (2-6) 
Changing variables x = £\/Ai in the integral (|2.ip gives 



p BL (v,t + At) = [ 

JQ 



PBL(£,t) 



exp 



4ira(£,VAt, t) 
(1 - exp 



ri-£,-a(£VAt,t)VAt 



4cr(CVAt,t) 



4a(£VAt,t) 



(2.7) 



The boundary layer solution has an asymptotic expansion in powers of V At 

p BL (v,t) ~ pgifa.i) + VAlp^lM + Mp%l(ri,t) + .... (2.8) 

Expanding all functions in (|2.7p in powers of V At and equating similar orders, we obtain 
integral equations that the asymptotic terms of (|2.8|) must satisfy. The leading order 0(1) 
term gives the Wiener-Hopf-typc equation on the half line 



P ( Bl(V,t) = 

for 77 > 0. The kernel 



y/4n<T(0,t) 









[ (v + 2 ' 




jexp 


4cr(0,t)_ 


+ exp 


4(7(0, t) _ 


} 





[ (v-0 2 ~ 




[ (v + H) 2 ] 


exp 


4ct(0, t) _ 


+ exp 


4cr(0, t) 



(2.9) 



(2.10) 



is an even function of 77 and £, i.e. K(r),£) = K{— 77, £) = K(r), — £) = A"(— 77, — £). Therefore, 
we extend Pg^i^t) to the entire line as an even function t) = P73i(~?ji)) 5 an d 

rewrite 



30 pgyy 

00 A/47T(7(0,t) 



exp 



4a(0,t) 



(2.11) 



for —00 < 77 < 00. The only solution of the integral equation (|2 . is the constant function, 
that is, PB L (f],t) = f(t), independent of 77. This follows immediately from the Fourier 
transform of (|2.1ip . whose right hand side is a convolution. 

Away from the boundary layer the solution admits an outer solution expansion 



Pour(y,t) ~p [ o UT {y,t) 



(2.12) 



where Pout satisfies the Fokker-Planck equation (| 1 .3(1 and the initial condition (|1.5p . Indeed, 
the integrals in (|2.1j) are of Laplace type with the small parameter At. For interior points 
y 3> V At the second integral, which represents only boundary interactions, is negligible 
relative to the first one. We change variables in (|2.ip by setting 



V 



y — x — a(x, t)At 
y /2a(x 1 t)At ' 



and extend integration over the entire line in the first integral and expand all functions in 
powers of V At. The resulting integrals are moments of the normal distribution. We obtain 

PAt(y,t + At) -p At (y,t) d[a(y,t)pAt(y,t)} d 2 [a(y,t)p A t(y,t)} 



At 



dy 



+ 



dy 2 



+ 0{VAt). 



The leading term in the expansion of pAt{y,t) is pQ UT (y,t), which therefore satisfies the 
Fokker-Planck equation (|1.3p . The initial condition (| 1 . 5|) is recovered from the Gaussian 
integral as At — > 0. The boundary condition that PouTiv^) satisfies can be determined 
only after the boundary layer is resolved by matching. The leading order matching condition 
of the boundary layer and the outer solutions is 



ton PsiM) =PouA°> t )- 



JO) 



Therefore 



p { £l(v,t)=pZ T M 



(0) 



(2.13) 



The matching condition at order yAt gives 
yPoliA ^) , „(i) 



dy 



+ PouA > t ) ~ PslX 7 ?.*) for »7 



which means that ?>sx,(?7, i) is asymptotically a linear function of rj, therefore the limit of its 
derivative is a constant. Thus the matching condition reduces to 



lim 

Tj — >00 



dp { B ] L (v,t) _ dp^ UT (0,t) 



drj dy 
The first order boundary layer term satisfies the integral equation 



exp 



P 



(v + 2 





[ + 




+ exp 




} 



(2.14) 



(2.15) 



o V4ttct(0, t) 



exp 



4cr(0,t) 



2a(0,t)7 ^4^(0, t) 





[ 




[ fa + a " 




jexp 


4cr(0,t)_ 


+ exp 


4cr(°)*). 


} 



4cr(0,t) 2 7 ^4^(7(0,*) 

4<r(0, t) 7 ^4^(0, t) 



£ U»7 - O 2 exp 



(»?-0 2 



(r? - f ) exp 



4a-(0,t) 



+ (jj + O 2 exp 



4cr(0,t) 



- (77 + f ) exp 



4cr(0,t) 
('7 + 2 



4cr(0,t) 



Evaluating explicitly the last four integrals in Q2.15|) and using (|2.13[) . gives 

d£- (2.16) 



V 4?rcr ( ' V 





[ (^-O 2 " 








jexp 


4a(0,i)_ 


+ exp 


4<r(0,t) _ 


} 



p 



»? 



^^ T (0^)erfcf 2V __ 



<r » (0 ' t) -° (0 '* ) .pgi T (0,t)«p 



VMM 

Differentiating (|2.16[) with respect to 77 and integrating by parts, we obtain 



4a(0,t) 



dr\ 
P 



2y/na(0,t) 
Setting 

we rewrite (|2.17|) as 

g(v,t) = 



y/4n<j(0, 1) Jo 
Poc/r(°'*) ex P 



#77 









[ (v + 2 ' 




jexp 


4a(0,f)_ 


— exp 


4cr(0,t) _ 


} 



4cr(0,t) 



q,(0,t)-a(0,t) (0) 
2^a(0,i) 3/2 p °£/TlU,t)r7exp 



4er(0,i) 



drj 



P 



Pc?i/ T (0,t)exp 



2^(7(0,0 



4o-(0,t) 



.(2.17) 



(2.18) 



(2.19) 



^4^(0, t) 7o 



g(£,t) < exp 



(r?-0 2 
4cr(0,7j) 



exp 



4cr(0,i) 



where 



P 



v/87rcr(0,t) 



Poc/r(M)exp 



8cr(0,t) 



q„(0,t)-ffl(0,t) (Q) fQ .s 



erf 



VMM 



(2.20) 



4cr(0,t) 



Since 4>(j],t) is an odd function of 77, we can define g(rj,t) for negative values as an odd 
function by setting g(r], t) = — g{— 77, rj) for r\ < 0. Then (|2.19[) can be rewritten as 



V47rcr(0, i) J- 



g(5,r;)exp 



4<r(0,t) 



which in Fourier space is 



g(k,t) 



<KM) 



l-exp[-cr(0,/j)fc 2 ]' 
Using the Wiener-Hopf method, we decompose 

g(k,t) = g+(k,t) +g-(k,t), 



(2.21) 



(2.22) 



(2.23) 



where g+(rj) = g{fl)x[o,co){v)i 3-0?) = g(v)X(-oo,o](v)- The Fourier transform g(k,t) exists 
in the sense of distributions, and g±(k,t) are analytic in the upper and lower halves of the 
complex plane, respectively. Taylor's expansion of 4>{k,t) in eq. (|2.20p gives 



0(k,t)=2ip { ° ) UT (O,t) 



PVV(M) 



[a y (0,t) - a(0,t)] }k + 0(k 3 ) as k -> 0. (2.24) 



The the non-zero poles of (|2.22[) split evenly between g + (k,t) and g-(k,t), and using 
g+(k,t) = —g—(—k,t), the pole at the origin gives 



.MM) = iPouA ^) 



p 



fj y (0,i) - a(0,t) 1 1 
cr(0, t) I k 



0{k) as fc->0. (2.25) 



Inverting the Fourier transform g + (k,t), by closing the contour of integration around the 
lower half plane, we obtain 



hm — ^ =Pou t(°'*) 

The matching condition (|2.14[) implies 



dy 



= PoutM 



P a y (0,t) - a(0,t) 

y/™(0,t) <r(0,t) 

a y (0,t) - a(0,t)\ 



P 



<r(0,t) 



(2.26) 



(2.27) 



Multiplying by cr(0, t) and rearranging, we obtain the radiation boundary condition 
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J(0,t) 



a(0,t)p^ UT (0,t)} - a(Q,t)p^ UT (0,t) = ™ ^ T (0,t). (2.28) 



9y 

Since p(y,t) — Pout^V^)i the reactive "constant" in (11.61) is 



«(*) = 



p^ioTt) 



(2.29) 



3. Numerical simulations in one dimension. The explicit analytical solution of 
the FPE (|1.3[) with the initial condition (| 1 . 5|) and the radiation boundary condition (|1.6[) 
for the case of vanishing drift (a = 0) and constant diffusion coefficient (a(x,t) — a) was 
first given by Bryan in 1891 [36] (see [33 §14.2, p.358]) 



p(x,t\ x Q ) 



1 



exp 



(a; - x ) 2 
Aat 



exp 



k { nix + x + nt) . 

- — exp < } erf c 

a a 



(x + x Q ) 2 
Aat 

x + xq + 2nt 



(3.1) 



'Aat 



The first term in (|3.1[) is the fundamental solution of (|1.3[) and (|1.5[) with a reflecting 
boundary condition, whereas the second term may be transformed into 



\J ira 3 t 



exp 



— r" exp 
a 



(X + Xq + Q 2 

4at 



which represents the density due to a line of exponentially decreasing sinks extending from 
— Xq to — oo. The method of Laplace transforming (|1.3[) with respect to t was later employed 
PP, [35] to obtain explicit analytical solution for the FPE (|1.3[) - (|1.5[1 with constants diffusion 
coefficient and (not necessarily vanishing) drift term a{x, t) = a 



p(x,t\x ) 



1 



\J\-KOt 

2n + a 



exp 



2(7 



■ exp 



(x — xq — at) 2 
M ^ 1 " X1) 

ax + k[x + x + (k + a)t 



axo (x + xq — at) 
a 



erfc 



4at 

x + x + (2k + a)t 



(3.2) 



Setting k = in (|3.2p reduces to Smoluchowski's [39] explicit analytical solution for a 
reflecting boundary with a constant drift term, while setting a = reduces to Bryan's 
solution (|3.1[) . 

We conducted several numerical experiments in which n = 10 7 trajectories were simu- 
lated according to the Euler scheme Q1.7P with the boundary behavior (|1.10[) . The diffusion 
coefficient was constant a — 1 and the reactive constant was k = 1, giving P — y^ir in 
cq. (|2.29p . The trajectories were initially located at xo = 1, their statistics were collected 
at time t = 1, and compared to the predicted p(x, t = 1 1 xq = 1). The convergence of the 
scheme was tested by using four different time steps At = 10 _1 , 10~ 2 , 10~ 3 , 1CT 4 . 

The first experiment corresponds to a vanishing drift a = 0. Figure 16.11 shows the 
convergence of the numerical scheme to the analytic solution (|3.ip . The rate of convergence 
of the numerical scheme to the analytic solution is V At. This is demonstrated, for example, 
by the survival probability 



Psur (X0,t) 



p{x, t I Xq) dx 



of finding the trajectory inside the domain at time t, that is, the probability that the 
trajectory was not absorbed prior to t. Integrating (|3.ip gives p sur (l, 1) = 0.77095 ... for 
a = k = 1. The survival probability is estimated numerically by the ratio of the number of 
survived (unabsorbed) trajectories n sur and the total number of simulated trajectories n = 
10 7 . Table l6Tl shows that the convergence rate of the estimated survival probability to its 
analytic value is \f~Ki as predicted by our boundary layer analysis. The statistical estimation 
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(variance) error due to the finite number of simulated trajectories is \Jp S ur{^ ~Psur)/ n — 
0.00013 . . ., which is an order of magnitude smaller than the smallest (bias) error obtained 
for At = 10" 4 (see Table |OJ. 

In the second experiment the drift term a = — 1 shifts the density leftwards, and causes 
more trajectories to react with the boundary. Figure 16.21 shows the convergence of the 
numerical scheme to the analytic solution (|3.2p . 

The final experiment corresponds to a reflecting boundary, P = k = and a constant 
non- vanishing drift towards the boundary a = — 1. We simulated n — 10 8 trajectories 
to obtain a finer resolution at the boundary. Figure 16.31 shows a comparison between the 
analytical solution (|3.2[) and the numerical densities for At = 10 _1 ,10 -2 . The no flux 
condition J = of a reflecting boundary together with (|1.4|) gives a negative boundary 
derivative p y (0, t) = —p(0, t) < 0. In particular, the analytic solution (|3.2[) satisfies p y (0, 1) = 
— p(0, 1) = —(2 + y/n) j (2-\/tt) ~ —1.06. The numerical densities, however, are flat at the 
boundary. Their first derivatives vanish at the boundary, as predicted in (|2.2|) and shown 
in Figure RP1 The first derivative changes from to 0(1) on an interval of length 0(y/At), 
manifesting a boundary layer behavior, though there is no such behavior in the density itself. 

4. Diffusion in M. d with partial oblique reflection at the boundary. We consider 
the d-dimensional stochastic dynamics 

x = a(x,t) + \/2B(t)w (4.1) 

in the half space 

Q = {x = (xi, x%, . . . , Xd) € K d : x\ > 0} 

where w is a vector of d independent Brownian motions and we assume that the diffusion 
tensor er(i) = B(t)B T (t) is uniformly positive definite for all t > s. The case of space- 
dependent diffusion involves many technically complicated calculations and will be consid- 
ered in a separate paper. We use henceforth the abbreviation <r(t) = cr. The radiation 
condition (|1.6[) becomes 

- J{V,t \x,s) ■ n = K(y,i)p(y,t \x,s), for y e dtt, x e £1, (4.2) 

where the components of the flux vector J{y, t\x,s) are defined by 

J k (y, t\x,s) = -[a k (y, t)p(y, 1 1 x, s)} + V A [a^ k p(y, t \ x, s)] , (4.3) 

where a^ ,k are the elements of the diffusion matrix cr. The Fokker-Plank equation for the 
pdf of x(t) can be written as 

dp(y,t\x,s) 



dt 



= -Vy • J(y,t\x,s) for all y,xetl. (4.4) 



x' = x + a(x, t)At + \/2B{t) Aw(t, At) fi, 



If x e f2, but 



the Euler scheme for (|4.1|) with oblique reflection in dQ. reflects the point x' obliquely in 
the constant direction of v to a point x" 6 il, as described below. First, we denote by x' B 
the normal projection of a point x' on 917, that ' — (x' ■ n)n. Then we write the 

Euler scheme for (|4.ip with partially reflecting boundary as 



x{t + At) = < 



for x 1 e Q 

w.p. l-P(x' B )\ r Al, if x'&Sl, (4.5) 



terminate trajectory w.p. P (x' B ) yAt, if x' £ fl. 
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The value of the termination probability P (x' B ) V At, that varies continuously in the bound- 
ary, is evaluated at the normal projection of the point x' on the boundary. The oblique 
reflection in the direction of the unit vector v (vi ^ 0) is defined by 



2x[ 

Vl 



(4.6) 



Note that x'{ = —x\ guarantees that the reflected point of a crossing trajectory is inside the 
domain fi. The fact that the normal components of x" and x' are of equal lengths makes 
the high-dimensional boundary layer analysis similar to that in one dimension. Normal 
reflection corresponds to v = n = (1, 0, . . . , 0). 

We note that for a point y e f2, we can write Pr{x" = y} = Pr{x' = y'}, where 



, 2y' ■ n 
y = y v 



is the oblique reflection of y' (see fig. 



Given y, equation (|4.7p defines y' as 



v = y 2— v. 

Vl 

As in the one-dimensional case, the forward Kolmogorov equation is 



(4.7) 



(4.8) 



p At (y,t + At) = [ — 

J Xl >0 (47T 



PAt{x,t) 



At) d / 2 Vdcter 
(l-P(y' B )7At)exp 



exp 



B{x + a(x,t)At,y) 



4At 

B(x + a(x,t)At,y')' 



+ 



4Ai 



dx, 



where 



B(x,y) = {x-y) T (T 1 (x-y) 



(4.9) 



(4.10) 



We construct a boundary layer of width O(yAt) in the normal direction to the boundary. 
The layer extends infinitely in the d — 1 directions tangent to the boundary 



Pbl{vi,V2, ■•-,Vd,t) = PAt{rfr"jA~b,y2, ■ ■ -,Vd,i). 

In other words, p B L(m n + Vbi *) = PAtiviVAtn + y B ,t), where y B = (0,2/2,2/3, 
As in the one-dimensional case, we assume the asymptotic expansion 

PBL(vin + y B ,t) ^ p ^l( m n + y B ,t) + y/Atp { g ) L {ri 1 n + y B ,t) + .... 

and substitute 



x = y B 



At£ 



in the integral (|4.9|) . We obtain 



p B L(vi n + y B , t + At ) = / 



>o 



exp 



£(£ + a(y B ,t)VAt,77in) 



exp 



1. 



B [£ + a(y B ,t)VAt, rjm 



PjnMin + y B + VAt£ B ,t) 
(47r) d / 2 % /deto : 

+ (l-P(y B )VAt)x 

Vl 



(4.11) 

■ -,yd)- 

(4.12) 
(4.13) 

(4.14) 



d£ + 0(At). 



ii 



We calculate separately the integral of the first and second terms in the braces. Substituting 
in the first integral 



transforms the domain of integration to 



(4.15) 



(4.16) 



where h = ^1/2^ is a unit vector, and a n = n 7 an = ||er 1//2 n|| 2 . Similarly, we transform 
the second integral by substituting z' = cr^ 1 /' 2 ( £ — mnH — — v ). Using the expansion 

V «i J 

(|4.12p . we obtain at the leading order the integral equation 



P { BL(mn + y B ,t) 
1 



(47r) d A j~„ . 
1 



Pbl ((Vi + V°n z ■ n) n + Vb, *) ex P 
Pbl i(-Vi + V°n z' -n)n + y B ,t)exp 



dz 



(47r) d A j z ,. n> jn n 

Integrating in the d — 1 directions orthogonal to n, yields 

1 f°° 

PBL(vm + y B ,t) = —= I Pbl((Vi + y/^v,)n + y B ,t)exp 
1 /"°° 

-7= / Psi {(-Vi + V^u)n + y B ,t)exp 

V47T Jjn 



dz'. 



21 



T 



1 



\/47rcr n Jo 



2>si (un + y B ,t) i exp 





(u-^i) 2 " 




(u + 771) 2 " 




jexp 


4cr„ 


+ exp 


4cr„ 


} 



This is the same leading order integral equation as that of the one-dimensional case (|2.9 
so the solution is independent of 771, and matching to the outer solution gives 



PBL^Vin + y B ,t)= p { ° } UT (y B ,t). 
To evaluate the 0{yfKt) terms, we expand in the first integral in (|4.14p 
B(£ + a(y B , t) VAi, 77m) = (£ - 77m) • cr" 1 (£ - 77m) + 



and in the second integral 



(4.17) 



(4.18) 



# ( £ + a(y B ,t)VAt, 77172, - — v 



At2a(y B ,i)-cr- 1 ^-7; 1 n). 



4 — 771 n H v ) ■ cr It — 771 n 1> 



M2a{y B ,t)-<r- 1 - m n^-v ) (4.1!)) 
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The 0(y At) contribution of the drift term for the first exponential term is 



.1 f PouAVByt) 

4 7 6>0 (47r) d / 2 VdeT^ 

(0) 

f 
4~ 



exp 



1 Pout (»B.*) n / ^ -1/2- f°° 
- = 2a(y B ,t)-o- n 

4 v47T J-vi/V^ 



we- tl2/4 dw = 



1 Pout{Vb^) , ,n f — »7l 

2 \/7rcr„ 4cr„ 



(4.20) 



The second exponential has the same contribution, so the overall contribution of the drift 
to the 0{\TKt) term is 



-a(y B ,t) ■ nexp 



/7T(T„ 



4ov, 



(4.21) 



(0) 
Pbl 



Now, we expand 

+ V^ z ■ ™0 n + Vb + V^* (ct 1/2 ^)b,*) {{ r h + V^ z ■ ™0 n + VB,t) 

'AtVpf L {{ m + ■ n)n + y B ,t) ■ {tr x l 2 z) B + O(M). 

Together with ()4.17|) . the expansion (|4.22|) reduces to 

Pbl + V^ z -n)n + y B + VA1 {ct 1/2 z) b 
P [ out (Vb, t) + VAt\7p^ UT (y B ,t) ■ {n 1 ' 2 z) B + 0{At). 



(4.22) 



Integrating as above, we obtain the 0{y At) integral equation as 
PBL(Vin + y B ,t) = 



V^On Jo 

p {v'b)Pout {Vb^) 



Pbl (v>n + y B ,t) < exp 





(u~T]l) 2 ' 




{u + Vi) 2 ' 




jexp 


4ct„ 


+ exp 


4ct„ 


} 



du 



exp 



{u + mf 



4ov, 



du 





" u 2 ' 


) exp 

' B 


4 _ 



du 



^/7T(J n 4cr, 



Differentiating with respect to 771 and integrating by parts (as was done in the one dimen- 
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sional case), we arrive at the integral equation 



r°° dpf L {un + y B ,t) 



exp 



V47tct„ Jo dn 

p (y'B)Pou T (yB^) 



^Pout(Vb^) ■{ — 



exp 



(u - rjif 



4a„ 



— exp 



(u + ?7i) s 
4cr„ 



du 



-»7i 



an t> 



771 exp 



4er„ 



crfc 



— v- 



2^/aZ 



B 



—. — a{y B ,t) ■ n- — exp 



-r,l 



4cr„ 



The Wiener-Hopf method requires the extension of the erfc function discontinuously as 
odd function, that is, to define erfc(x) = sgn(x) erfc(|x|). Following the calculations of the 
one dimensional case, it remains to determine the small k behavior of the Fourier transform 
of erfc(cc). Using 



as an 



crfc 



V 



2,/a 

we obtain, as in (|2.24[) . 

p (y'B)PouT (Vb^)^ 



exp{ikn} dn <~ 2ik j erfc ( — - — J r\dn — 2ika n , (4.23) 

2\/ (7n I 



4>(k) ~ 2ik< 



- 2a n Vp { o> UT {y Bl t) ■ 



an v 

a n 2vi 



PoUT{VB,t) a {yBi~t) n \ a s k^O. 



Therefore, 



Um d P ( B L { m n + y B ,t) = 
dn 



Tji — *00 



P(v'b)Pqut (VB,t) 



2Vp^( Vfl ,t) 



an v 

On 2v\ 



(o) , ,,a(y B ,t)-n 
' Pout l« b > l ) 



Combining with the matching condition 

Hm dp { BL {nm + y B ,t) = dp ( ^ UT {y B ,t) 

we obtain 



n — >oo 



dn 



dn 



(4.24) 



dn 

P(yB)PouT (VbJ) 



hxa r . 



^p { ouT(y B ,t) 



an v 

<j n 2vi 
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The requirement that the pdf of the limiting diffusion process satisfies the Robin boundary 
condition leads to the only possible choice 



crn 



(4.25) 



Otherwise, we obtain an oblique derivative boundary condition. Since y' B — > y B as At — > 0, 
we obtain the Robin boundary condition 

-JouT(y B ,t) ■ n = ^PouAvb^) ■ (rn-p { o}j T {y B ,t)a{y B ,t) ■ n 
_ P{Vb) Pout (Vb^) 



The reflection direction v of crossing trajectories is the co-normal direction crn. Normal 
reflection (i.e., replacing v by n) gives rise to the boundary normal flux if and only if n is 
an eigenvector of the diffusion tensor cr. The limit of the outer solution as At — > is the 
solution of the Fokker-Planck equation (|4.4p with the radiation boundary condition 

-J(y,t)-n = K (y)p(y,t) for y E <9f2, (4.26) 

where the reactive "constant" is 

n{y) = — . 4.27) 

Note that normal reflection will not recover the normal flux of the radiation condition if n 
is not an eigenvector of cr. 

5. Numerical simulations in two dimensions. To illustrate the co-normal reflec- 
tion law (|4.25p in the Euler scheme (|4. 5[) - (14. 7|) in the half plane x > 0, we ran several 
numerical experiments. The simulations show the convergence of the pdf of the numerical 
solution to that of the FPE with the radiation boundary condition (|4.26[) - (|4.27[) . Unlike the 
one-dimensional case, no explicit solution of the anisotropic Robin problem for the FPE in 
the half plane is available, so we compare the statistics of the simulated trajectories with 
a numerical solution of the FPE. The latter is constructed by the stable Crank-Nicolson 
scheme on lattice points, where in each time step the sparse linear system is solved by the 
conjugate gradient method. 

In all numerical experiments the initial point is (%o,yo) = (0.3,0) and the statistics is 
collected at time T — 0.5. We choose the reactive constant n = 1 and the diffusion matrix 
B in gUD 



B 



0.3 0.4 
1 



which gives the anisotropic diffusion tensor 

, T /0.25 0.4 



cr = BB 



0.4 1 



We simulate n = 10 7 trajectories with time steps At = 10 _1 , 10~ 2 , 10 -3 , 10" 4 in each 
experiment. 

In the first experiment the drift vanishes, a — 0. The normal n = (1,0) and the 
co-normal crn = (0.25, 0.4) point in different directions. The simulated trajectories are 
reflected in the co-normal direction according to the prescription (|4.25[) . The simulated and 
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the numerical solution of the FPE give the marginal densities shown in Figures 16.51 and 16.61 
Figure [ITU shows the marginal density of x(T), 



p(x,T\x ,y ) = p(x,y,T\x ,yo)dy, 

J — OO 

while Figure RTTTTI shows the marginal density of y(T), 

p(y,T\x 0: y ) = p(x,y,T\x ,y )dx. 



Table 16.21 gives the computed survival probability and indicates the convergence rate. 

We illustrate the importance of using the correct reflection law in the second experiment, 
in which the simulated trajectories are reflected in the normal direction n = (1, 0). Clearly, 
the marginal density of x{T) coincides with that of the first experiment, because both 
oblique and normal reflections have the same ^-coordinate (see (|4.6|) 1. However, the plot of 
the marginal density of y{T) differs significantly from that in the previous experiment. It 
is apparent from the comparison to the numerical solution of the FPE that the simulation 
does not recover the Robin boundary condition in the limit At — > (see Figure |6~7)) . Note 
that the peak of the density is at y > 0, though the reflection is normal. This is due to the 
anisotropy of the diffusion tensor, which causes the probability flux density vector to have 
a positive y component. 

In the third experiment the drift is the constant vector a = (—1, 0) and the diffusion ten- 
sor is as in the first experiment. The density is shifted toward the boundary (see Figure [6781 
and Figure [679]) . The results are summarized in Table [6~3l 



6. Summary and Discussion. We have defined a diffusion process with partially 
reflecting boundary as a limit of Markovian jump processes generated by the Euler scheme 
for the dynamics in a half space with partial absorption of exiting trajectories and partial 
oblique reflection in the boundary. We derived an expression for the radiation constant in 
the Robin boundary condition for the one-dimensional Fokker-Planck equation for the case 
of diffusion with variable drift and diffusion coefficients, as a function of the absorption 
probability. We found that the Euler scheme for a diffusion in a half space with variable 
drift and constant anisotropic diffusion has to be reflected in a particular oblique direction 
in order to recover the Robin boundary condition. Also for this case we found the radiation 
" constant" as a function of the local absorption probability on the boundary. We found a 
boundary layer of width 0{\fKt) in the pdf of the Euler scheme and solved the boundary 
layer equation, which is of Wiener-Hopf type. 

The boundary layer of PAt(y, t) makes the calculation of the boundary flux non-trivial. 
The net boundary flux of the simulation profile PAt{y,t) is 



, . I Py/Ai f° , [°° , . f (x-yf 

JAt{0,t) = lim — / dy I pAt{x,t)exp< - 

At^o At ^/AnaAt J-oc Jo { 



AnaAt 



dx, (6.1) 



which is the probability of the trajectories that propagate per unit time out of the domain, 
discounted by the probability of trajectories returned into the domain by the partially re- 
flecting Euler scheme. Changing the order of integration and then changing the variable of 
integration into z — x/2\/aAt gives 

- Ja*(0, t) = ( erfc(z)p At (2zV^At, t) dz = (0, t) + O(VaI). (6.2) 

Jo V 7 !" 

This straightforward calculation of the flux gives the correct radiation constant, provided 
that 

P { bIM=P { outM- (6-3) 
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The latter, however, depends on the mode of reflecting a trajectory from x' outside to x" 
inside the domain. We have shown that for x" = —x' the provision holds, however, for other 
schemes, e.g., x" — —ax' (a 7^ 1), the provision (|6.3p fails in general, though (|6.2p still 
holds. On the other hand, the differential form of the flux, (ll.4p . has to be obtained from 
(|6.ip in the limit At — > 0, which is not the case for PAt{y,t), though it is for pouT(y,t). 
This shows up in spades in the multi-dimensional case, because although (|6.3p holds for any 
direction of reflection, yet the differential form of the flux is obtained in the limit only if the 
correct direction of oblique reflection is chosen. 

The generalization of the multi-dimensional case to domains with curved boundaries 
and to a variable diffusion tensor tr(x,t) is not straightforward and will be done separately. 
Note that if the diffusion tensor is constant, but un-isotropic, a local orthogonal mapping of 
the boundary to a plane converts the diffusion tensor from constant to variable, as can be 
seen from Ito's formula. However, as mentioned in Section [TJ in the most common case of 
constant isotropic diffusion, our result extends to domains with curved boundaries, because 
the mapping leaves the Laplacian unchanged, though the drift changes according to Ito's 
formula. In this case the vector v coincides with the normal n. 
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At 


Ti sur 


Psur W> sur / Tl 




7253450 


0.0456 


io- 2 


7577156 


0.0132 


10~ 3 


7670969 


0.0039 


10~ 4 


7698523 


0.0011 



Table 6.1 



Survival probability: the difference between the analytic value of the survival probability p S ur = 
0.77095... and its numerical estimation n sur /n decreases by roughly VlO whenever At is decreased by 
an order of magnitude. (Parameters: a = K = xo = t = l, a = 0, n = 10 7 j 



At 


Ti sur 


Psur flgur/fl 


IO" 1 


5986662 


0.0814708 


10- 2 


6449991 


0.0351379 


10~ 3 


6707318 


0.0094052 


10- 4 


6775672 


0.0025698 



Table 6.2 



Survival probability for a = 0. The third column lists the error between the numerical value of the 
survival probability p SU r = 0.6799545 from the solution of the FPE and its estimate n sur /n from the 
simulation. The error decreases by about \/T0 whenever At is decreased by an order of magnitude, indicating 
the convergence rate %/At of the simulation. 



At 


Ti sur 


Psur W> sur / Tl 


io- 1 


2541947 


0.1180946 


10- 2 


3399528 


0.0323365 


10~ 3 


3632622 


0.0090271 


10~ 4 


3693905 


0.0028988 



Table 6.3 

Survival probability for a = (—1,0). The third column lists the error between the numerical value of 
the survival probability p sur = 0.3722893 from the solution of the FPE and its estimate n aur /n from the 
simulation. 
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Fig. 6.1. No drift: The analytical solution f3.1\) (Magenta), and the three numerical densities At = 
10~ 1 (Blue), At = 10 -2 (Green), At = 10 -3 (Red) approaching it from below. The numerical density of 
At = 10 — 4 is not shown, because it is difficult to distinguish it from the analytic density. (Parameters: 
a = n = XQ = t = l, a = 0, P = y/n, n = 10 7 ) 
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Fig. 6.2. Drift, a 
1CT 1 (Blue), At = 10" 
xq = t = 1, P = y/n, n 



-1: The analytical solution (3. 21) (Magenta) and the numerical densities At 



(Green), At 
W 7 ) 



10" 



(Red) that approach it from below. (Parameters: 




Fig. 6.3. Drift, a = —1, reflecting boundary P = k = 0: The analytic solution f3. 2]) (Red) and the 
numerical densities At = 10" 1 (Blue) and At = 10 — 2 (Green) with n = 10 s simulated trajectories to obtain 
a finer boundary resolution. (Parameters: o = n = XQ = t = l) 
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Fig. 6.4. A simulated trajectory can get from x to y in a single time step At in two different ways: 
(i) directly from x to y, without crossing the boundary, and (ii) by crossing the boundary from x to y' and 
reflection in the oblique direction v with probability 1 — P(y' B )y/ At to y. The reflection law (^.5[ )- f4~?l ) 
satisfies y' x = —yi . 
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Fig. 6.5. The marginal density of x(T) with no drift and correct oblique reflection (the first experiment). 
The numerical solution of the FPE (blue) with grid size Ax = 0.01 and estimates from the simulation of 
n = 10 7 trajectories with time steps At = 10 _1 , 10 -2 , 10 -3 , 10 -4 . 
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Fig. 6.6. The marginal density ofy(T) with no drift and correct oblique reflection (the first experiment). 
The numerical solution of the FPE (blue) with grid size Ax = 0.01 and estimates from the simulation of 
n = 10 7 trajectories with time steps At = 10~ 1 , 10~ 2 , 10 -3 , 10 -4 . 
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Fig. 6.7. The marginal density ofy(T) with no drift and with normal reflection (the second experiment). 
The numerical solution of the FPE (blue) with grid size Ax = 0.01 and estimates from the simulation of 
n = 10 7 trajectories with time steps At = 10~ 1 , 10~ 2 , 10 -3 , 10 — 4 . 
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At = 0.1 
At = 0.01 
At = 0.001 
At = 0.0001 

Numerical Ax = Ay = 0.01 




Fig. 6.8. The marginal density of x(T) with drift a = (—1,0) and correct oblique reflection (the third 
experiment). The numerical solution of the FPE (blue) with grid size Ax = 0.01 and estimates from the 



simulation of n = 10 trajectories with time steps At = 10 ,10 ,10 



, 10" 
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Fig. 6.9. The third experiment (a = [— 1,0] T , correct oblique reflection): y-marginal densities. The 
numerical solution (blue) is compared to four simulated ones (with time steps At = 10 _1 , 1CT 2 , 1CT 3 , 10~ 4 ). 
n = 10 7 . Resolution: Ax = 0.01. 
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